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ABSTRACT 

We investigate the secular dynamics of a planetary system composed of the parent star and 
two massive planets in mutually inclined orbits. The dynamics are investigated in wide ranges 
of semi-major axes ratios (0.1-0.667), and planetary masses ratios (0.25-2) as well as in the 
whole permitted ranges of the energy and total angular momentum. The secular model is 
constructed by semi-analytic averaging of the three-body system. We focus on equilibria of 
the secular Hamiltonian (periodic solutions of the full system), and we analyze their stability. 
We attempt to classify families of these solutions in terms of the angular momentum integral. 
We identified new equilibria, yet unknown in the literature. Our results are general and may be 
applied to a wide class of three-body systems, including configurations with a star and brown 
dwarfs and sub- stellar objects. We also describe some technical aspects of the semi-numerical 
averaging. The HD 12661 planetary system is investigated as an example configuration. 
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1 INTRODUCTION 



days, about thirty extrasolar multi-planet systems have been 
eqj Many of them seem either locked in or close to low- 



Nowada 
detecte( 

order mean motion resonances (MMRs). Moreover, there is a class 
of the so called hierarchical systems jLee & Peale 2003 ) which can 
be characterized by relatively small ratio of semi-major axes. Their 
planetary orbits are well separated and far from collision zones, 
hence the long-term, qualitative dynamics of such systems may be 
effectively investigated with secular theories. The Hamiltonian of 
a hierarchical system can be averaged out over mean longitudes 
which play the role of fast angles. In the regime of small eccen- 
tricities and inclinations, this approach leads to the well known, 
classical Laplace-Lagrange (L-L) secular theory (Murray & Der-| 
|mott 2000 ). It relies on the expansions of the disturbing function in 
power series with respect to eccentricities and inclinations which 
are small parameters of the problem. However, many multi-planet 
hierarchical systems do not satisfy the assumption of small eccen- 
tricities, and the L-L theory may fail. 

Still, to deal with the observed diversity of orbital configura- 
tions, the secular theories relying on high-order expansions of the 
perturbations are used, e.g., the series in eccentri city (e.g., | Mur 



|^& Dermott 2000, Rodriguez & Gallardo 2005[ [Libert & Hen- 
Irard 2006 2007a lb 1 1 Vera s & Armitage 2007) or exnansion to the 



Tard 2006 20 07a|b| |Veras & Arm itage 2007 ) or expansion to the 
third order in the ratio of semi-major axes a, known as the octu- 
ple theory pord et al.|2000l|Lee & Peale|2003] . This theory can 



* E-mail: c.niigaszewski@astri.uni.torun.pl 
t E-mail: k.gozdziewski@astri.uni.torun.pl 

^ See Jean Schneider's Extrasolar Planets Encyclopedia http://exo planet.eu| 
for frequent updates on the discoveries and orbital parameters 



be also generalized to higher orders ( Migaszewski & Gozdziewski| 
|2008| and references therein). The analytical expansions are partic- 
ularly suitable for studies of hierarchical systems. Moreover, they 
are usually valid only in limited ranges of the orbital parameters, 
and special cases (like resonant configurations) must be treated 
individually. The alternative, recently developed quasi- analytical 
theory relies on averaging the perturbing Hamiltonian numerically 
(Michtchenko & Malhotra 2004 , Michtchenko et al.|2006j . In this 
work, we are heavily inspired by these papers and their idea of 
the semi-analytical technique. Because the method does not re- 
quire any expansion of the perturbing Hamiltonian, basically, it has 
no limitations inherent in the analytical theory. For instance, with 
a help of this technique, Michtchenko & Malhotra, (2004) found 
new, non-classic feature of the secular dynamics of coplanar sys- 
tem of two planets (the so called non-linear secular resonance in 
the regime of large eccentricity). In the later work, [Michtchenko | 
[et al. I f 2006) consider more general three-dimensional (3-D) secu- 
lar model of two-planet system, and present a systematic approach 
helpful to investigate the global dynamics of such configurations. 

As an example to study, we choose the HD 12661 planetary 
system ( [Fischer et aIl|200T| [2003| [Butler et al.][2006l ). The dis- 
covery paper ( [Fischer et al.[|2001| ) announces two Jovian planets 
on well separated orbits with semi-major axes of ^0.8 au and 
~ 2.8 au, respectively, and of moderate eccentricities. We analyzed 
the most recent, publicly available data from the catalogue of [But-[ 
[ler et al.[ ( |2"006 ) and ( Wright et al.|2008j ), using the A^-body model 
of the radial velocities (RV) and the so called hybrid minimization 
( [Gozdziewski & Migaszew ski'2006V The results of our analysis of 
the RV data published in ([Butler et al. 2006 ) are illustrated in Fig.[T] 
The best fit solution yielding (Xv)^^^ ^ 108 and an rms ~ 7.5 m/s 
is marked with a crossed circle in the dynamical map in terms of 
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the Spectral Number (SN) ( [Michtchenko & Ferraz-Mello|[200T] ). 
The SN is the fast indicator making it possible to distinguish be- 
tween chaotic and regular planetary configurations. The osculating 
elements of the best fit solution at the epoch of the first observation 
are given in caption to Fig.[T] In this figure, we mark the semi-major 
axis and eccentricity of the outer planet derived from an ensemble 
of fits within la of the best fit solution. Clearly, the available data 
already constrain orbital elements of the outer planet very well. The 
dynamical maps reveal orbits well separated from the low-order 
MMRs. Two most prominent MMRs within the vicinity of the best 
fit are 19:3 and 13:2 MMRs, respectively. Moreover, the best- fits 
within la confidence level span the region of small eccentricities 
in which the resonances are very narrow. 

Hence, the HD 12661 system fits well assumptions of the sec- 
ular theory. This system has been studied already in a few papers: 
with the direct numerical integrations (J i et al. 2003 ), with the ana- 
lytical octupole theory of hierarchical systems ( Lee & P eale 2003 ), 
with mapping of the phase space by fast indicators ( [Gozdziewski| 
[2003 a| ), and with the classic analytical theory that relies on ex- 
pansions of the perturbation in eccentricity ( Rodriguez & Gallardo] 
|2005[|Liberr& Henrard 2006 ) . All the cited works assume that the 
HD 12661 system is coplanar and oriented edge-on. However, we 
should keep in mind that a major limitation of the Doppler tech- 
nique lies in the ambiguity of orbital inclinations, which cannot 
be well determined by far. The observational windows are still rel- 
atively narrow, and to remove the inclination degeneracy, several 
orbital periods of the outermost orbit are required. Moreover, the 
recent formation theories do not fully predict mutual inclinations 
in multi-planet systems. We cannot be certain yet whether the com- 
mon assumption of coplanar orbits really holds true. Likely, many 
different forming scenarios are possible. For instance, the migra- 
tion in low-order MMRs may end up with systems characterized by 
large mutual inclinations ( |Thommes & Lissauer|200 3 ). The dynam- 
ical relaxation of initially dense planetary systems of giant planets 
( [Adams & Laughlin|2003| ) may lead to scattering events which pro- 
duce wide distribution of the mutual inclinations. Indeed, recent 
simulations of [Veras & Ford| ( |2008| > revealed that the outer planet 
of the HD 12661 system undergoes large oscillations for nearly all 
of the allowed two-planet orbital solutions. These authors conclude 
that it might be the effect of a perturbation of planet c, perhaps due 
to strong scattering of an additional planet that was subsequently 
accreted onto the star. Moreover, we stress that the inclination of 
the HD 12661 system is still unknown, hence the understanding of 
basic features of its 3-D dynamics seems also important. This in- 
triguing system is an excellent candidate for tests and numerical 
experiments regarding non-coplanar configurations. Moreover, be- 
cause we attempt to study the secular 3-D dynamics globally, our 
results are general and valid for much wider class of three-body 
systems. The secular theory which we consider here, covers plane- 
tary systems with different masses and semi-major axes ratios, and 
the full range of mutual inclinations. 

The plan of this work is the following. In Sect. [2] we re- 
call the general mathematical model of the 3-D two-planet sys- 
tem. Subsection 12.11 is devoted to a short technical overview of 
the averaging approach and some computational details that may 
be useful in a practical implementation of the method. To make 
the paper self-contained, we also recall the notion of representa- 
tive planes, and energy levels calculated for fixed values of the 
total angular momentum integral (Sect. |3]). To illustrate the pre- 
cision of the semi-analytic approach, we compute the Poincare 
cross sections, and demonstrate chaotic behavior of the secular sys- 
tem ( ^Michtchenko et al.|2006j ). The main results are described in 



Sect. |4] which is devoted to the analysis of the existence and bi- 
furcations of equilibria in the secular, spatial problem of two mu- 
tually interacting planets. In particular, we detect and investigate 
closely a few families of these equilibria in a wide range of plan- 
etary mass ratio, fi G {0.25,0.5, 1,2}, and the semi-major axes ra- 
tio, a e {0.1,0.2,0.333,0.667}. The results are general and valid 
as long as the partition of the Hamiltonian onto the Keplerian, 
integrable term, and the small perturbation is reasonable. In that 
part, our work extends the paper of [Libert & Henrard] p007b] ). 
After introducing non-singular canonical variables, they investi- 
gate the existence, stability and bifurcations of stationary solutions 
emerging from the equilibrium at zero-eccentricities, the so called 



Lidov-Kozai resonance (|Lidov|1961 Kozai|1962 


Lidov & Ziglin 


19761 [Thomas & Morbidem|1996|[Innanen et al. 


19971 [Kinoshita 


& Nakai[|2007( to mention a few papers in an endless list of ref- 



erences) which was found and intensively investigated in the re- 
stricted three body problem. The full three-body problem in the 
Hill case (hierarchical configurations) was also intensively stud- 



ied by many authors (e.g., Krasinsk y|1972 [T974[[Lidov & Ziglin 
[1976, Ferrer & Osacar|1994f [Miller & Hamilton|2002| [Fabry cky 
|& Tremaine 2007|. These works rely mostly on the second order 



expansion of the secular Hamiltonian in the semi-major axes ratio 
(the quadrupole approximation). In the present work, we focus on 
two aspects of the problem: 

• we consider the unrestricted problem in wide ranges of semi- 
major axes ratio a, up to 0.667, and mass ratio ^, 

• we study equilibria of the full secular Hamiltonian; the semi- 
analytic averaging helps us to compute the secular perturbations 
beyond convergence limits of the usual power series expansions. 

Thanks to the quasi- analytic averaging, we found new families of 
equilibria of the secular 3D planetary problem which unlikely may 
be detected with the help of perturbation techniques. We also study 
the Lyapunov stability of these solutions in detail (or to an extent 
permitted by technical limits of the semi-analytic algorithm). 



2 THE 3-D DYNAMICS OF TWO-PLANET SYSTEM 

The Hamiltonian of the three-body planetary system, expressed 
with respect to canonical Poincare variables (see [Laskar & Robutel] 
1995 , Michtchenko et al. 2006 ) has the form of: 




Keplerian part 



k^mim2 ^ Pi P2 

A ^ mo ' 
direct part indirect part 



(1) 



where r/ denote the position vectors relative to the star, p/ - their 
conjugate momenta relative to the bary center of the full three-body 
system, A = |ri — is for the distance between planets, mo - is 
the mass of the parent star; mi , m2 - are for the masses of the plan- 
ets (also index / = 1 is for the inner planet, and / = 2 for the outer 
planet). We denote also the mass parameters fJl = (mo + m/) and 
the reduced masses through = (1/m/ + 1 /mo)~^ . Under the as- 
sumption of mi <C mo (or, more generally, for small enough per- 
turbations of Keplerian orbits), the Hamiltonian of the system, 
is a sum of the Keplerian term (which would be integrable in the 
absence of mutual interactions between planets) and the interaction 
term with the so called direct and indirect terms. 

Alternatively, the dynamical state of the system, (r/,p/) may 
be represented through the mass-weighted canonical angle-action 
variables of Delaunay ( [Murray & Dermott^2000j : 
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Figure 1. The dynamical map of the edge-on, coplanar HD 12661 system in the (flc,^c) -plane, for the best-fit solution to the RV data published in |Butler| 
|et a l.'20()6^>. Large values of the Spectral Number (SN) marked in yellow indicate continuous spectrum of the fundamental frequencies of the system and 
strongly chaotic motions, small SN (black) means discrete frequencies and ordered motions. Positions of low-order MMRs are labeled. The best-fit solution, 
yielding (Xv)^^^ ^ 1-08 and an rms ~ 7.46 ms~^ in terms of parameter tuples (m [mj], a [au], e, CO [deg], M [deg]) including planetary mass, semi-major 
axis, eccentricity, the argument of pericenter and the mean anomaly at the epoch of the first observation, ro=JD2,450,83 1.608380, is the following (2.34, 0.831, 
0.361, 296.24, 126.86) for planet b, and (1.84, 2.888, 0.021, 52.66, 66.18) for planet c, respectively. The original errors are rescaled by adding stellar jitter of 
3.5 m/s in quadrature. The mass of the parent star is 1.11 Mq. Solutions within la level of the best fit are marked with yellow circles, fits with marginally 
worse (Xv)^^^ are marked in red. 



(2) 



Mi — the mean anomaly, L/ = y^/L^~ai, 

CO/ — the argument of pericenter, G/ = Li ^\—e^ 
Q.i — the longitude of node. Hi = Gi cos //, 

where a/ denote semi-major axes, ei - eccentricities, and // stand 
for inclinations; {Li,Gi,Hi) are the conjugate momenta. The trans- 
formation between (r/,p/) and the set of Delaunay eleme nts 
{ai,eiJi,Q.i,(Oi, Mi) may be found in ( [Ferraz-Mello et al.|2005| ) or 
( |Morbidelli|2002| ). If orbits are far from strong MMRs and collision 
zones then Hamiltonian in Eq.[T]can be averaged out with respect 
to the mean anomalies playing the role of fast angles, and then we 
obtain the secular Hamiltonian i^ec which approximates the long- 
term, slow variations of the mean elements. 

To make the paper self-contained, we recall basic facts on the 
secular 3-D problem of two planets. We follow Mich tchenko et aL\ 
(|2006 ) and Libert & Henrardj(|2007b|. The averaged i^^c does not 
depend on Mi,M2, therefore the conjugate actions (Li,L2) are 
constant. After the partial reduction of nodes, i^ec depends on AO. 
only, and not on Q^i and separately. After the canonical trans- 
formation ( [Michtchenko et al.,2006j : 

(coi,Gi) (0)1, Gi) 

(C02,G2) (C02,G2) 

(ei = i {ai + a2),Ji=Hi +H2) 



(3) 



{a2,H2) (02 = \ {^1-^2). J2=Hi -H2) . 

The secular Hamiltonian does not depend on 9i, therefore Ji = 
\C\= const, where C is the total angular momentum of the system. 
Moreover, 62 = 7i/2 = const (in the Laplace frame, 202 = AO. = 
±71, after Jacobi's elimination of the nodes), and: 

J2 = {Gl-Gl)/Ji. 

For fixed angular momentum Ji and secular energy i^ec^ the av- 
eraged system can be reduced to two degrees of freedom. In- 
stead of Ji , we may use the so called Angular Momentum Deficit, 
AMD = Li + L2 — /i . The AMD is a measure of the system non- 
linearity ( |Laskar|2000| ). Coplanar and circular orbits have the min- 
imum of AMD = 0. In configurations with large AMD, crossing 
orbits are possible and they are very unstable. 

Because the secular Hamiltonian still depends on many pa- 
rameters, the global analysis of the long-term dynamics are com- 



plex. To simplify the study of their basic properties, we fix par- 
ticular values of integrals and/or orbital elements. We choose the 
semi-major axes and masses as the primary parameters of the sec- 
ular model. Then Li and L2 are their (scaled) representation. The 
maximum of AMD is equal to Li +L2, hence we introduce the nor- 
malized Angular Momentum Deficit: 

which is an uniform and non-dimensional representation of AMD. 
Relative to the Laplace plane, Cx = Cy = 0, = C, hence we have: 



Li ^ 1 — ^^cos/i +L2 ^1—^2 cos/2 = Jli 
Li \—e\ sin/i —L2 \j\—e\ sin/2 = 0. 

Also the mutual inclination of orbits /^ut = h-^h- Thus, cos /^ 
cos/i cos/2 + sin/i sin/2 cos All, or, alternatively, 

T2_r2_ r2 

C0S/mut(^l,^2) = 



2G1G2 



C0S/i(^l,^2) 



COS/2(^l,^2) 



/f + G 



2/1 Gi 



/f + G^ 



(4) 



(5) 



(6) 



2/1 G2 

Because C = Q > 0, the above relations are singular for Ii = l2 = 
7i/2or^l=^2 = l (when Gi = G2 = 0). A boundary of the mani- 
fold of permitted motions for a given Ji=C (or AMD), semi-major 
axes and planetary masses ratio, can be defined through /^ut = 0, n. 
It can be also shown that when AMD is fixed and /i 2 < 7i/2 then 
the mutual inclination at the origin {ei = 0, ^2 = 0) is maximal. We 
will denote such value by /q from hereafter. 

The dynamics of the secular system are expressed through so- 
lutions to the following canonical equations of motion: 

dco/ ai^^ec dOi d^, 



dGi 



dt 



1,2, 



(7) 



where (c0i,C02) are canonical angles and (Gi,G2) are canonical 
momenta. Having only the numerical approximation of i^ec (see 
below), we must solve Eqs. [7] numerically. For that purpose, we 



© 0000 RAS, MNRAS 000, 000-000 



4 C. Migaszewski and K. Gozdziewski 



may use a suitable integrator, for instance, the 4-th order Runge- 
Kutta scheme ( [Press et al.|1992| |. The partial derivatives appearing 
in the right-hand side of the equations of motion, are calculated 
with the mid-point rule ( [Press et al.|1992| ). Moreover, to calculate 
precisely the Hessian of i?^ec which is required to determine the 
stability (see Sect. [3.2| ), we are forced to use higher order approxi- 
mations of the second order partial derivatives. 

2.1 The semi-analytical averaging 

The problem is now to average out the Hamiltonian, Eq. [T] We 
calculate: 

^^^^=(2SF h Jo (8) 

where 9{ is the Hamiltonian of the problem expressed through the 
canonical Delaunay elements. For small enough perturbations, the 
Keplerian part of (i^) depends on constant Li only and does not 
affect the secular evolution of the system. It can be shown that the 
average of the indirect part of Hamiltonian equals to a constant 
( jBrouwer & Clemence|1961| ). In the non-resonant case, we have to 
average out the direct part of disturbing Hamiltonian only. 

The analytical calculation of apparently trivial integral, Eq. [8] 
is in fact a difficult problem. Usually, the Hamiltonian is expanded 
in power series with respect to appropriate small parameter (eccen- 
tricity, inclination or semi-major axes ratio). Then with the help 
of a suitable canonical transformation, we can "remove" particu- 
lar terms of the Hamiltonian. However, the secular series converge 
for relatively small values of parameters. Instead, as we mentioned 
akeady, the secular Hamiltonian, Eq. [8] can be computed numer- 
ically, without troublesome power series expansions. This bright 
idea of [Michtchenko & Malhotra ( 2004 ) is quite simple to apply. 

Apparently, to compute integral in Eq.js] we must evaluate 
in a discrete grid of the mean anomalies. That would imply multi- 
ple (and in fact unnecessary) solution of the Kepler equation. To get 
rid of this problem, we can change the variables under the double 
integral using the well known expressions relating the mean {Mi), 
true {fi) and eccentric CE^) anomalies, respectively. To express the 
double integral through the true anomalies, we differentiate the Ke- 
pler equation Mi = "Ei — Ci sin 1^, with respect to Mi^^, and then 
we find that dMi = Ji dfi, where: 

Ji = MeiJi) = {l-efy^^ {l^eiCosfi)-\ i=l,2- (9) 
The secular Hamiltonian has the following form: 

^sec--^^ / / 7dhdh. ^ = ^Jlh- (10) 
(ItzY Jo Jo 

We may also express the double integral through eccentric anoma- 
lies that leads to even simpler expressions for functions Ji. Next, to 
calculate the integral in Eq. [TO] we apply an adaptive-grid integra- 
tion algorithm that relies on the Gauss-Legendre quadrature of the 
64-th order. The adaptive algorithm is forced by large variability of 
the integrand function. To illustrate that issue, we analyse a few typ- 
ical examples shown in Fig. [2] The left-hand contour plots in this 
figure are for the shape of direct term of ^ (Eq. [!} multiplied by 
JiJ2^ in the (/i 5/2) -plane. These plots are computed for differ- 
ent values of eccentricities and AG5 = 05 1 — 052, where 05 1^2 are the 
longitudes of periastron. In this experiment, the system is copla- 
nar. In the top-left panel of Fig. [2] (see its left-hand part), which 
corresponds to relatively small eccentricities, is weakly varying 
function of (/i 5/2). But for large eccentricities, it may have narrow 



extrema in some parts of the (/i,/2) -plane (see the bottom-right 
contour plot in Fig.|2]). In these areas, to reach a desired accuracy, 
the integral must be computed on a dense grid of the (/i , /2) -plane. 
However, in other parts of the grid, such a large number of quadra- 
ture nodes is not necessary and, under the requirement of fixed ac- 
curacy, it would cause significant CPU overhead. Thus, the optimal 
computation of the double integral is possible with the non-uniform 
grid in the (/i,/2) -plane, following an idea of adaptive quadra- 
tures [see, for instance, [Press et aL] ( [1992| )]. In the right-hand part 
of Fig. |2] we illustrate the steps of our simple adaptive mesh in- 
tegration by appropriate divisions of the integration subintervals. 
Typically, the number of divisions is small but in some parts of the 
(/i 5/2) -plane, it may be as large as 8-9, in order to obtain the rela- 
tive error of 8 ~ 10~^^ in two subsequent steps of the integration. 



3 EQUILIBRIA IN THE 3-D SECULAR PROBLEM 

According to the classic methodology of Poincare, to understand 
the dynamics, one should investigate whole families of solutions. 
Isolated orbits in the phase space tell us little on the global proper- 
ties of the system. The most simple class of solutions that can be 
investigated efficiently in any two degree of freedom Hamiltonian 
system are equilibria defined through algebraic equations: 



Typically, one tries to find the phase-space coordinates of these 
equilibria, their number and bifurcations as well as to determine 
their Lyapunov stability (at least, the linear stability). The analy- 
sis of the existence and bifurcations of equilibria in the secular 3D 
system are quite complex because they depend on many parame- 
ters (AMD, the total energy, particular orbital elements, masses of 
planetary companions). Hence, to investigate such solutions glob- 
ally, we have to choose a proper representation of the phase space 
regarding these parameters. Moreover, to avoid limitations of the 
analytical approach, the whole analysis should be done numeri- 
cally, by the semi-analytical averaging. Hence, a reduction of the 
dimension of the phase space is critically important. 

3.1 The representative planes of the energy 

To simplify the search for equilibria of i^ec» we choose a specific 
two-dimensional plane of initial conditions that makes it possible 
to represent the stationary solutions in the 4-D phase-space of the 
secular system. We follow [Michtchenko et al.[ p006| ) and [Libert| 
|& Henrard| ( [2007bt . The representative plane of initial conditions 
(the 5-plane from hereafter) should have common points with each 
phase trajectory of the secular system. In |Michtchenko et al.[p006| l, 
the 5-plane is defined through: 

= {ei cosA05 X ^2cos2c0i}, 

where ei^2 ^ [O5 1] and angles (A05, 2coi) are fixed to pairs of angles 
(0,0), (71,0), (0,7r), and (7C,7c), respectively. In that notion, the 5- 
plane comprises of four subsets of points which coordinates span 
the range of ei cosA05 G [—1, 1] and ^2cos2c0i G [—1, 1]. 

Subsequent panels of Fig. [3] show generic views of the ^Pm- 
plane derived for different values of integral and the same pri- 
mary parameters, (a,]u). In particular, these plots are drawn for i^^ec 
levels which are found numerically as solutions to ^sec — Eo = 0, 
where Eo is a fixed value, for the following values of the semi-major 
axes and masses ratios: a = ai/a2 = 0.333, ju = mi/m2 = 0.5, 
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Figure 2. The left-hand panels are for contour levels of function (iTd in the (/i , /2)-plane, computed for the coplanar two-planet system and orbital parameters: 
mo = 1 Mq, mi = 1 mj, m2 = 3 mj, <3i = 1 au, ^2 = 3 au. Eccentricities and A03 are different at each panel: the top left-hand panel is for ei = 0.1, ^2 = 0-2, 
Acn = 0, the top right-hand panel is for e\ = 0.6, 62 = 0.5, A03 = 7r/2, the bottom left-hand panel is for e\ = 0.4, 62 = 0.5, A03 = 71, the bottom right-hand panel 
is for e\ = 0.6, 62 = 0.7, ACH = Ti/2. Panels in the right-hand column illustrate the AMR-like division of the integration domain, as depending on the variability 
of the integrand function. 



-0.85 




Ci cosAtD 



e, cosAtD 



cosACD 



Figure 3. Generic plots of i^ec in the ^Pm -plane of (^i cos A03, 62 cos2coi ) obtained for a = 0.333, fji = 0.5 and = 0. 10 (a) and ^ = 0.25 (b) and Si = 0.85 
(c), respectively. Small Roman numbers label families of stationary solutions identified in this work (see the text for more details). Solid thick line marks the 
collision line of orbits defined through ai{\ ^6\) = a2{l — 62). Shaded regions mark mutual inclinations. Panel (a) [0°,50°] (whit6), [50°, 60°] (light gray), 
and ^ 60° (dark gray), respectively. The limit values of inclinations for panel (b) are, 0°, 95° and 105°, respectively; for panel (c): 0°, 140° and 155°. 



and = 0.10 (the left-hand panel), = 0.25 (the middle panel), 
and = 0.85 (the right-hand panel), respectively. According with 
the general construction of the 5-plane, it is divided by four quad- 
rants and, for a reference, labeled with Roman numbers in the left- 
hand panel of Fig. [3] quadrant I (AG5 = 0, COi =0), quadrant II 
(AG5 = 71, COi = 0), quadrant III (AG5 = 0, COi = 7i/2), and quad- 
rant IV (AG5 = 71, COi = 7i/2). 

Alternatively, we also use another definition of the 5-plane: 

= [d sincoi X ^2sinC02 : C0i,C02 = ±7i/2;^i^2 ^ [0, 1)}, (12) 

(see [Libert & Henrard|2007b| . The ^p9-plane helps to avoid a dis- 
continuity of the levels of 5£ec at the x-axis. In fact, that plane 
carries out the same information as the negative (y < 0) part of 
the ^pM-plane. Obviously, pairs of angles of the ^Ps representa- 
tion: (COi = +7l/2,C02 = -7l/2), (COi = +7l/2,C02 = +7l/2), (COi 

71/2. 0)2 = — 7i/2), (coi = — 7i/2,C02 = +7i/2), correspond to the 
following pairs of angles in the ^Pm representation: (AG5 = 0, 2cOi = 
71), (AG5 = 71, 2cOi = 71), (AG5 = 71, 2cOi = -71 ), (AG5 = 0, 2cOi = -7l). 
Hence, two bottom quadrants of the fP^ -plane are equivalent to 



quadrants IV and III of the -plane. Two upper quadrants of the 
^p5-plane are their central reflections with respect to the origin. It 
follows from the definition of coordinate axes through ^/sin() and 
^jCOs() functions (where ij =1,2). Apparently, the fP^-plane con- 
tains redundant information. However, the energy levels are con- 
tinuous in this plane and their interpretation is easier than in the 
^pM-plane [see also ( Libert & Henrard 2007b )]. The central projec- 
tions of quadrants III and IV can be obtained by reversing signs of 
COl and CO2 (or measuring angles in opposite direction). 

We define one more 5-plane, which makes it possible to obtain 
a smooth representation of quadrants II and I of the ^pM -plane: 

= {^iCOSCOi X ^2COSC02 : C0i,C02 = 0,7i;^i^2 ^ [0, 1)}- (13) 

Because we are interested in possibly global and transparent repre- 
sentation of the equilibria in the secular problem (see below), we 
will use not only the primary notion of the 5-plane by |Michtchenko| 
|et al.| (2006 ) but also the two other definitions. 

An important observation which is very helpful to justify the 
choice of the 5-planes for the search for equilibria, is the symmetry 



© 0000 RAS, MNRAS 000, 000-000 



6 C. Migaszewski and K. Gozdziewski 



of i^ec with respect to the characteristic plane. It can be shown as 
iefin( 



follows. For the defined above pairs (co^ CO2) of the 5-plane: 



d^ I(co?,co0) 



:0. 



(14) 



:co?,co?) 



Indeed, from the general formulae of the secular Hamiltonian ex- 
pressed by Fourier series we have: 



h,l,m (ai , ^2 , ^1 , ^2 , A , ^2 ) cos <I), 



k,Lmi 



k,l,me{ 



where kj,m are integers, are coefficients of the expansion, 
and ^k,l,m = kodi + /CO2 + mAQ. is the generic angle argument 
of the expansion. (Further, we shall assume that the series con- 
verge). According with the analytic properties of the Fourier expan- 
sion, indices k and / must have the same parity ( (Brumberg|1995[ 
[Michtchenko et al.1|2006| . Also, after the Jacobi's elimination of 
nodes, AO, = ±7r. Now, the derivatives of ^Hgec over co/ (Eq. 14 1 are: 



dG/ 
dt 



hk,l,m (^^ 1 , ^^2 , ^ 1 , ^2 , A , ^2 ) sin <I), 



LLr 



'k,l,m 



k,l,me{ 



and because coefficients can be considered as functionally 

independent, the derivatives may vanish only when all sin<I>yr. = 
0. This is only possible when ^k,l,m = n E Z, hence, when 
^coi + /CO2 = ±{n — m)7i, for any integers kj of the same parity, 
and when co/ = ±7i/2,0,7i. That also means, that 2coi = 0,7i and 
AG5 = G5i — 052 = 0,^- 

The zeros of the derivatives of the secular Hamiltonian over 
CO/ may be also deduced geometrically, relying on the symmetry 
of interacting mean orbits. The mean orbits may be understood as 
material elliptic rings (the Gauss approximation), which interact 
gravitationally. The potential of interaction has symmetries with 
respect to the particular angles AG5,2coi or (C0i,C02) which define 
the 5-plane. Points (G9, GS) in the 5-plane, fulfilling conditions: 



= 0, 



(G?,G?,co?,co?) 



(15) 



may be identified with stationary solutions of the secular problem. 
We solve the above equations with respect to unknown (GpG2) 
or, (^p ^2) VSiirs of fixed angles (cOp CO2) in the given quadrants 
of the 5-plane and for fixed C. Hence, the notion of the 5-plane is 
particularly suitable for the analysis of equilibria. 

Figure[3]reveals numerous stationary solutions labeled accord- 
ingly with the quadrant of the -plane and a letter labeling a spe- 
cific type (a family) of solutions. The equilibria appear as local ex- 
trema (or rather as elliptic or quasi-elliptic points) or saddle points 
of i^ec in the 5-plane. At these critical points, the derivatives with 
respect to all phase variables must be equal to zero. After fixing the 
(a, /i)-pair, ei and 62 may be varied in ranges permitted by constant 
C = Ji (or AMD). The thick curve is for the boundary of the energy 
level defined for a given value of AMD. The eccentricities and mu- 
tual inclination are coupled again through Ji (or AMD). To indicate 
boundaries of the mutual inclination permitted for a given range of 
{€1,62), the regions in which the mutual inclination is grater than a 
prescribed value are shaded. We mark a few such shaded regions in 
the 5-plane (lighter shade means smaller mutual inclination). The 
mutual inclinations at their boundaries are quoted in the caption to 
Fig.[3](also in captions to other plots of the 5-plane). 

To avoid the geometric singularity of the equations of motion 
at the origin of the 5-plane and at the coordinate axes (x = ei =0, 
y = 62 = 0), WQ follow ijLibert & Henrard| (j2007bf , and introduce 
the following non- singular, canonical variables: 




-0.45 0.45 
ej cosACiJ 



■0.5 0.5 
ej cosACiJ 



Figure 4. Levels of the secular Hamiltonian in the quadrupole approxima- 
tion (top panels), and in the octupole approximation (bottom panels). These 
plots are obtained for a = 0.333, jn = 0.5 and = 0.1 (the left-hand pan- 
els), and = 0.25 (the right hand panels). Compare with the semi-analytic 
model in Fig. [3] The shaded areas mark parameters for which the expansion 
of fH^Qc in terms of a would diverge. 



Pi = ^/2{Li-Gi) cos CO/, qi = a/2(L/-G/) 



smcOr. 



1,2. 



We denote x = {pi,qi,P2iQ2) from hereafter. These non-singular 
variables are convenient for a quasi-global continuation of station- 
ary solutions in the 5 -plane. 

Finally, to show the relevance of the semi- analytic averag- 
ing, we calculated the energy levels in the 5-plane when only the 
quadrupole (~ a^) and octuple (~ a^) terms of the perturbing 
Hamiltonian are accounted for. These terms are averaged analyt- 
ically. The results are illustrated in four panels of Fig. [4] which are 
derived for the same values of a = 0.333 and ju = 0.5 as in Fig. [3] 
Panels in the top row are for the quadrupole-order secular theory, 
panels in the bottom row are for the octupole theory. The left-hand 
plots are for = 0.1, the right-hand plots are for = 0.25. Shaded 
areas mark regions of the parameter plane which lie beyond the 
limit of convergence of the expansion of i^ec in oc, and obviously 
we cannot obtain there a proper representation of equilibria solu- 
tions. The quadrupole term leads to exactly symmetric view of the 
5 -plane — in fact, the quadrupole Hamiltonian does not depend 
on AG5. The octupole approximation fits much better to the semi- 
analytic secular model (compare with Fig. [3^,b), nevertheless the 
energy levels are still significantly distorted and some features are 
missing at all; for instance, there is no quasi-elliptic point over the 
collision line in quadrant II (see Fig.[3]D); instead, we may found a 
false saddle solution close to the border in quadrant I. Although the 
tested configuration has relatively small a = 0.333, in such a case 
both analytic approximations of i^ec introduce artifacts which can 
be only avoided by an application of the semi-analytic averaging. 



3.2 Lyapunov stability and critical inclinations 

The stability of equilibria may be investigated with the help of 
Lyapunov theorem (see, e.g., Markeev "1978; "Khalil 2001). If the 
Hamiltonian is positive (negative) definite function in a neighbor- 
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hood of an equilibrium xq, then the equiUbrium is Lyapunov sta- 
ble. At a stable equilibrium, the parameters of the averaged sys- 
tem are constant, hence the orbital elements do not change in 
the secular-time scale, and the orbital configuration evolves in the 
short- time scale only. In the 3-D secular problem, that is equiv- 
alent to conditions for local extrema of i^ec in the phase space. 
We recall that such extrema must appear as elliptic points in 2-D 
plots of the 5-plane. We should also remember that these quasi- 
elliptic points may be in fact related to saddles in two remaining 
and "hidden" dimensions of the phase space. To determine, whether 
the secular Hamiltonian is sign definite function of the phase vari- 
ables in the neighborhood of a critical point, we compute its Hes- 
sian, H2 = ^^i^ec/^x^ at the equilibrium, and then we determine 
whether it is sign-definite matrix. 

As an illustration, we show two plots of the sub-determinants 
of Hessian IHI2 in Fig. |5] which are computed at the origin as func- 
tions of the mutual inclination of circular orbits, /q- The top-panel 
is for a = 0.333, id = 2 and the bottom panel is for a = 0.667, 
jj = 2. For some particular values of io, which are related to the 
given value of J^, the sub-determinants may vanish and then we 
cannot determine sign-definiteness of the Hessian. For instance, if 
a = 0.333 then it takes place for /q close to - 45°, 70°, 130°, and 
155°, respectively. 

In fact, these values of /q are related to bifurcational incli- 
nations of the "trivial" equilibrium at the origin ( Krasinsk y(1972| 
1 1974) and changes of its stability and global topology of i^ec- 
The later work gives explicitly their values in terms of parame- 
ter P = L1/L2 ~ id^/a which were calculated for the second or- 
der secular Hamiltonian (the quadrupole term). For a reference, 
the vertical lines in Fig. [5] mark the bifurcations derived with the 
quasi- analytic theory (thin lines), and with the quadrupole Hamil- 
tonian (thick, dashed lines). Bifurcatio nal values of io ar e labeled 

the "-h" 



with ^1 2 3 4- Following terminology of Krasinsky 



1974 



sign means that the bifurcation of the ongm leads to nontrivial so- 
lution of the positive type, the "-" sign means nontrivial solution 
of the negative type. The positive type solutions are characterized 
by ^1,2 = 0?^. hence bifurcations take place in the ^Pc-plane; the 
negative type equilibria (001,2 = ±7i/2) appear after bifurcations in 
the fP^-plane. An inspection of Fig.[5]reveals, that the bifurcational 
inclinations may be very different in both theories, and it may be 
particularly well seen for a = 0.667 (bottom panel of Fig. [5]). In the 
later case, the bifurcational values of inclination are clearly split- 
ted, and the bifurcations takes place for different Si (or Ji =C). 
Note that they only depend on P in the quadrupole theory, and on 
fd and a separately in the full model. Actually, angles = and 



2 are degenerated in the quadrupole theory (note that the oc- 



tupole theory breaks the symmetry). All that means that the topol- 
ogy of the phase space must be different in the two secular theories. 

When the Hamiltonian evaluated at a critical point is not a sign 
definite function then the analysis of stability become much more 
difficult than in the case of an extremum. In general, only the linear 
stability of the equilibrium can be determined relatively easy. We 
accomplish that by solving the eigenproblem of matrix A of the 
linearized equations of motion. The variational equations in terms 
of new canonical variables y, where x = xq + y : 



dy 
&t 



:Ay, A = m2(xo), 



and I is the symplectic unit. In general, for a conservative Hamilto- 
nian system, A has complex eigenvalues 
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Figure 5. The sub-determinants of the Hessian of the secular Hamiltonian, 
deti H2, det2 IHI2, deta H2 and det4 H2 evaluated at the origin of the 5-plane. 
The sub-determinants are plotted with gradually shaded curves as functions 
of Si or Iq, starting with the black curve for deti H2, and lightest gray curve 
for det4El2. The sub-determinants are expressed in relative units. The ver- 
tical lines mark the bifurcational mutual inclinations of orbits at the origin 
(io) that correspond to det/ EI2 = for one (or more) i= 1,2,3,4. These in- 
clinations are labeled with I^^^ a- Thin, dashed, vertical lines labeled at the 
top are for the quadrupole order theory, and the thin solid lines labeled at 
the bottom are for the semi-analytic theory. See the text for more details. 



We can find them easily as the roots of symmetric characteristic 
polynomial p(X) = det(A — XE) = 0, where E is the unit matrix. 
It is well known that the necessary and sufficient condition for the 
linear stability is fulfilled if Xi = ±ia/ are purely imaginary and 
matrix A is diagonalizable (a/ are the characteristic frequencies). 

In the case of two-degree of freedom conservative Hamilto- 
nian systems, we can apply the theorem of Arnold-Moser (e.g., 
[Meyer & Schmidt|1986| ) to conclude that equilibria which are lin- 
early stable are generically Lyapunov stable. However, there is no 
such implication if the characteristic frequencies are involved in 
resonances up to the 4-th order, i.e., when pci + qa2 = for < 
\p\ + kl ^ 4, with p,q eZ, or when coefficients of the Birkhoff's 
normal form of the Hamiltonian expanded near the equilibrium 
fulfill a particular condition involving a/ [see ( [Meyer & Schmidt] 
[1986| ) or ( |Markeev,1978j for details]. In resonant cases, we should 
examine each particular normal form of the polynomial expansion 
of the Hamiltonian with respect to variations y. This can be done 
with the help of constructive theorems by Markeev and Sokol- 
skii [see, e.g., jMarkeevj ( [1978] ); jSokolskiij ( |1975| ) or |Gozdziewski[ 
( ,2003b] ) for an example application of these theorems, and refer- 
ences therein]. Moreover, because high-order expansions are re- 
quired, such an extensive study is hardly possible because we must 
average out i^ec (^nd calculate its derivatives numerically. A pre- 
cise enough determination of the second order derivatives becomes 
very difficult. Hence, we are forced to limit the stability analysis to 
the linear, non-resonant case. Nevertheless, recalling the implica- 
tions of the Arnold-Moser theorem, a study of the linear stability 
provides valuable information on the generic Lyapunov stability. 
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3.3 A general view of the 5-plane 

While in Sect. 4 we describe the results regarding new families 
of stationary solutions found in this paper in some systematic way, 
here we refer to generic properties of the 5-plane identified through 
many numerical experiments. Fixing J^, we can obtain typical 
views of the 5-plane which are shown in three panels of Fig.|3] The 
le ft panel of Fig. [3] illu strates configurations recently investigated 
in |Michtchenko et al.| (i2006) and |Libert & Henrard| ( |2007b| ). We 
can see a maximum of i^ec at quadrant IV (coi = 7c/2, 0)2 = 7r/2) 
of the 5 -plane. It corresponds to equilibrium marked with IVa and 
known as the Lidov-Kozai resonance, with the analogy to the re- 
stricted problem (Lidov 1961 ; Kozai 1962 ). In the vicinity of equi- 
librium IVa of the non-restricted problem, angles coi and CO2 librate 
around 7c/2. Simultaneously, these librations of coi 2 are related to 
large-amplitude, anti-phase variations of the eccentricity of the in- 
ner orbit and of the mutual inclination. This mechanism may lead to 
strong instability. We observed it already in the case of hierarchical 
two-planet configurations ( [Gozdziewski & Konacki|2064| . 

Due to discontinuity of the fP^-plane at the x-axis, it is difficult 
to follow the evolution of geometric structure of the L-K resonance. 
Instead, the and ^Pc representative planes are more convenient 
for that purpose, particularly near the origin. A sequence of plots 
shown in Fig.|6] reproduces the analytical results of [Libert & Hen-| 
|rard| ( [2007b] ) which were obtained for a = 0.1 and fji = 0.25. For 
= 0.01 (the left-hand panel of Fig. [6|, the origin is stable, per- 
mitting mutual inclination of circular orbits /q ~ 30°. With increas- 
ing = 0.03, the inclination grows, and for /q ~ 43°, the stable 
stationary point become unstable and bifurcates. Three new solu- 
tions appear: one is unstable and two are stable. This phenomenon 
may be called the L-K bifurcation. At the bifurcation point, some 
sub-determinants of IHI2 are equal to zero and the stability cannot be 
determined (see Fig. [5] and the previous Section for details). With 
further increase of J^, the L-K resonance centers move toward large 
values of ei (see the third panel in Fig.|6]plotted for = 0.06) and 
approach ei ~ 1 for = 0.08 (see the l ast, fourth panel in Fig. [6] ). 
While we refer to the analytic work of ( [Libert & Henrard|2007b| , 
these authors did not follow the L-K equilibrium for this value of 
J^. The semi-analytical algorithm makes it possible to continue the 
family of L-K solutions up to such a value, for which we observe 
new bifurcations of the equilibria. From each bifurcation of stable 
L-K equilibrium emerge three new solutions: one linearly stable (a 
saddle point in the 5-plane) and two elliptic points. One of them is 
Lyapunov stable, the other one is unstable. The elliptic points ap- 
proach ei and moderate 62- The solution at the origin bifurcates 
the second time but it remains unstable (note that it appears as an 
elliptic point in the 5-plane) and two unstable equilibria (saddles of 
i?^ec) located at {ei ~ 0,^2 > 0) also appear. 

In the second plot of the 5-plane (see the middle panel of 
Fig. [3]), we consider a configuration with a = 0.333, ju = 0.5 for 
= 0.25. We can recognize the L-K resonance after a bifurcation: 
the bottom-left quadrant of the 5-plane reveals a local extremum 
labeled with IVa and a saddle point IVb-. In remaining three quad- 
rants, we can find also other new equilibria labeled with Ila, lib. 
Ilia, Illb, respectively. Curiously, the maximum marked with Ila 
lies beyond the geometrical crossing line of orbits defined implic- 
itly through ai(l±^i)=a2(l — ^2) - Finally, in the right-hand panel 
of Fig. [3] we draw the ^pM-plane for large = 0.85 which lead to 
a discontinuity of the energy plane in the regime of large ^2- In this 
case, the inclination reaches very large values. 

These examples indicate that the 3-D problem is much more 
complex and rich in dynamical phenomena than the co-planar prob- 



lem of two planets. We recall that in this case ( [Michtchenko &| 
Malhotra 2004 ), the phase space of the secular system is spanned 
by librations of AG5 around (mode I), librations of AG5 around 
71 (mode II), and circulations of AG3. There is also possible the 
so called non-linear secular resonance (the true secular resonance) 
which is present in the regime of moderate and large eccentricities. 

Generally, the equilibria are not isolated in the parameter 
space of (;W, a) and the AMD integral. Yet, as the example of the 
L-K resonance demonstrates, the stationary solution may evolve 
in the parameter space, they can bifurcate, and may change their 
stability. Hence they form families of solutions and their behavior 
depends on a complex way on problem parameters. To investigate 
these families, we require a continuation method for determining 
bifurcations of the equilibria and their stability. 



4 PARAMETRIC SURVEY OF EQUILIBRIA 

To resolve the families of equilibria, we apply a simple continua- 
tion method with respect to as the primary parameter. For fixed 
parameters a and /li, ^ [O5 !]• We increase this quantity by small 
steps, and we compute the secular energy map in the 5 -plane. An 
inspection of the characteristic plane makes it possible to detect the 
origin and development of basic dynamical structures. In particu- 
lar, we can determine critical values of for which new equilibria 
(represented by elliptic or saddle points in the 5-plane) appear (see, 
e.g., Fig.[6|. Having an overall view of the 5-plane, we may follow 
a given solution along some path in the parameters space with the 
help of a minimization algorithm (see below). 

Here, we show two example sets of the 5 -plane derived for 
two pairs of (a,^). Figures |7|{9] are for {a^jn) = (0.333,0.5), while 
figures 10 ■ 1 1 are for {oL,ju) = (0.2, 2). We start to look more closely 
at the first set of the energy diagrams. Figure|7]comprises of a num- 
ber of panels derived for varied values of J^. Each value of is re- 
lated to the mutual inclination at the origin, /q- Shaded ares in these 
plots mark ranges of the mutual inclination permitted by the given 
and fixed J^. Lighter shadings encode smaller mutual inclinations. 
Thin black curves encompass the region of permitted motion ac- 
cording with /mut = 0, n. We show the 5-plane defined as (Pm (see 
Fig. [7]) as well as the iP^-plane (Fig.jSj and fPc-plane (Fig. [9]). We 
can see very clearly that the fP^^c-planes provide a continuous rep- 
resentation of energy levels. 

A sequence of energy diagrams shown in Figs.|7j|9]helps us to 
understand the development of a few families of equilibria found 
in this work. We start with = 0.03 (the top left-hand panel). In 
this case, the origin is the global extremum (the maximum) of the 
secular energy. This corresponds to the well known classic zero- 
eccentricity equilibria investigated in detail in ( [Krasinsky|[1972[ 
P^74l [Libert & Henr^[2007b| ). When = 0.1 (the next panel 
in the top row), we see a saddle at the origin and a maximum in 
quadrant IV of fP^-plane, which can be better seen in the ^p5-plane. 
The extremum can be identified with the Lidov-Kozai resonance. 
Clearly, in that case, the neighboring trajectories characterized by 
librations of angle COi around 7i/2. We can also notice a non-classic 
feature, regarding the non-restricted model of the L-K resonance: 
close to the quasi-elliptic point, also angle AG5 may librate around 
71 (it means that CO2 librates around 7i/2). This effect is possible for 
compact systems. 

When = 0.18 (the next panel in the top row of Figs. [7]- 
|9]), new structures appear: a saddle at the origin of the ^Pc-plane 
(see appropriate panel in Fig. [9]) with two elliptic points close to 
the ^1 =0 axis, as well as an elliptic point above the collision line 
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Figure 6. Levels of the secular energy (J^ec) at the ^P^-plane defined through (ei sincoi, g2sinc02) with 0)1,0)2 = ±7t/2, for a = 0.1, ja = 0.25 and are 
(from the left-hand panel to right-hand panel): 0.01 , 0.03, 0.06, 0.08, respectively. Shaded areas indicate ranges of the mutual inclination: (0° , 20° , 30°) for the 
left-hand panel, (0° , 50° , 60° ) and (0° , 70° , 80° ) for the inner panels, respectively, and (0° , 95° , 1 05° ) for the right-hand panel. 



(marked with thin lines). The change of topology of i^ec is also vis- 
ible in the respective panels of the ^p5-plane. As we already noticed, 
it is related to the second bifurcation of the origin. Simultaneously, 
the center of the L-K resonance moves toward large ei. When the 
= 0.25 (the top right-hand panel of Fig.jSj we observe a further 
development of the structure around the origin and a bifurcation 
of the extremum identified with the L-K resonance onto a saddle 
point and two new elliptic points appearing in the regime of large 
ei. These structures are particularly well seen in fP^c-planes, re- 
spectively, as illustrated in Fig.|8]and Fig. [9] 

A similar analysis might be carried out for the second set of 
parameters (/i, a). We do not present it in detail, however an inspec- 
tion of the energy diagrams shown in Fig.[TO] and Fig. [TT] reveals 
that the sequence of plots ends in a different dynamical situation: 
even for large J^, the region of permitted motion remains closed. 
Obviously, the evolution of equilibria appearing for the second pair 
of parameters (a,^) is different from the first case. 

Finally, we consider one more experiment devoted to a com- 
parison of the results derived with the help of the octupole theory 
and the quasi- analytic averaging algorithm. Figure [12] illustrates a 
few families of equilibria derived for a = 0.333 and ju = 2. Black 
and red filled circles are for the semi-analytic theory, blue and green 
filled circles are for the octupole Hamiltonian approximation. The 
red and green points indicate equilibria found beyond the formal 
limit of convergence of the expansion of i^ec in oc. We may notice 
significant differences between some branches of equilibria already 
in the regime of moderate 62- There are also solutions permitted by 
the octuple theory, represented by green horizontal branches, which 
are absent in the quasi- analytic model. This test confirms that the 
study of equilibria in compact systems benefit from the application 
of the semi-analytic (basically exact) averaging. 

Hence, we should follow a more systematic procedure. Once 
we identify a solution of a given family for fixed pair of (a,/i), we 
may continue that family by searching for the zeros of the right 
hand sides of the equations of motion, Eqs. [15] This task may be 
accomplished by minimization of the norm of the partial derivatives 
of the secular Hamiltonian. To speed up the minimization, we apply 
the fast Levenberg-Marquardt algorithm ( [Press et al.|1992| ). Simul- 
taneously, we examine the stability of solutions which are found 
after the L-M algorithm converged. Perhaps a more elaborate algo- 
rithm of the continuation of the equilibria might be applied, never- 
theless, even with such a simple approach, we are able to identify 
a few families of solutions that are, to the best of our knowledge, 
unknown in the literature. Finally, Figs. []3[fT5] illustrate the results 
of the continuation globally. Each set of panels is derived for fixed 
(yW, a) chosen as combinations of parameters a = 0.2, 0.333, 0.667 



and ju = 0.25, 0.5, 1, 2, respectively. The continuation of the fam- 
ilies of equilibria is done in the whole possible range of J^. We 
plot /mut. ^1 and 62 of the found stationary solutions, as functions 
of (and, if the value of permit circular orbits, as a function 
of io). Columns in each group of diagrams are for particular quad- 
rants of the 5-plane. Note, that we skipped panels for quadrant I 
because in that quadrant we found only one family la of unstable 
solutions for a limited range of /i ^ 1 (see Sect. |4.2| for details). 
Simultaneously, Lyapunov stable (or linearly stable) equilibria are 
marked with large filled circles, and unstable solutions are marked 
with small filled circles. Families of equilibria are classified accord- 
ingly with the quadrant of 5 -plane in which they appear, and they 
are labeled with corresponding Roman numbers. Let us note that 
the red filled circles indicate equilibria found beyond the formal 
limit of convergence of the expansion of i^ec in a. 

In this way, we can obtain quite a deep insight into the secu- 
lar equilibria and their stability in wide ranges of the primary pa- 
rameters. Below, we describe the identified families of stationary 
solutions in more detail, and we try to characterize the associated 
dynamical behaviors of the secular system. A likely position of the 
HD 12661 system in the diagrams Figs. [T3}fT5] may be deduced 
from its currently known orbital elements, a ~ 0.3 and jn ^ I. For 
relatively small ~ 0.1, the system might be found in the top, 
right-hand panels of Fig. [14] in the regime of the L-K bifurcation, 
still with moderate ei ^0.3 and 62 ^0.1. 



4.1 Family at (ei , 62) = (0, 0) 

The stationary solution at the origin {ei = 0,62 = 0) (see an ex- 
ample in Fig. |6]) was investigated in detail by [Libert & Henrard] 
(2007b) for a = 0.1, A/ = 0.25, i.e., for relatively distant orbits (or 
typically hierarchical configuration). In our classification, this fam- 
ily is marked with "0" in all stability diagrams of Figs. T3]- 15] So- 
lution s of this family are also studied in detail by |Krasinsky |l972| 
1974| ) in terms of the quadrupole approximation of i^ec, and we al- 
ready did many references to these works and its results. Here, we 
start to follow more closely the evolution of family "0" for {a,]u) : 
(0.333,0.5) with respect to Jl. Similarly, Fig. 



10 



and Fig. 



11 



veal topology of the 5-plane and the evolution of zero-eccentricity 
equilibria for a different pair of parameters, {oi,]u) = (0.2, 2). 

When the mutual inclination remains relatively small (see 
Figs.[7}|9]), the zero-eccentricity equilibrium is Lyapunov stable be- 
cause it corresponds to the maximum of i^ec- When increases, 
a bifurcation of this solution appears for /q ~ 43° (as mentioned 
already, the L-K bifurcation). Inspecting the ^p5-plane (Fig. [5]), we 
may notice that at the bifurcation point a new family IVa appears. 
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Figure 7. Levels of the secular energy (i^ec) at the ^P^-plane defined through {ei cosACH, g2cos2c0i) with A03 = 0,71, 2cOi = 0,71, for a = 0.333, /u = 0.5 
and varied R (or the mutual inclination of circular orbits, io). Values of R are (counting plots from the left to the right and from the top to the bottom): 
0.03,0.1,0.18,0.25,0.33,0.4,0.46,0.50, respectively. Each panel has two shaded regions indicating ranges of the mutual inclinations (corresponding to the 
fixed above of J^): (20°, 30°), (50°, 60°), (70°, 80°), (95°, 105°), (120°, 130°), (140°, 150°), and (130°, 140°), (120°, 130°), respectively. 




and it remains Lyapunov stable up to extremely large ei. After ei 
reaches the limit of 1 (simultaneously, the mutual inclination is 
close to 7i/2), the family bifurcates again: the family IVa remains 
on a branch with large ei while a new branch of solutions (family 
IVb-) can be continued to large mutual inclinations with simultane- 
ous decrease of ^i. Family IVb- can be regarded as the retrograde 
case of the L-K resonance with /mut > 7c/2. Note that families IVa 
and IVb- are quasi- symmetric with respect to /mut ~ 7i/2, regarding 
the eccentricity and inclination of the inner orbit. This symmetry is 
more "exact" for smaller a (hierarchical configurations). We notice, 
that the quadrupole term approximation leads to exact symmetry of 
the equilibria (see Fig.|4]and the relevant comments in Sect. 4.3). 

We may also discover unstable family Illb, which has the el- 



liptic point located in quadrant III as well as a saddle of family 
IVb-i- corresponding to linearly stable solution. In the later case, 
the neighboring trajectories are characterized by librations of coi 
around n/2 and also (within a limited vicinity of the libration cen- 
ter), by librations of CO2 around — 7i/2. A further increase of leads 
to a shift of solutions IVa, IVb-, and Illb towards the border of per- 
mitted motions. Finally, for > 0.4, the zero-eccentricity equilib- 
ria vanish at all, and the energy plane is divided onto two distinct 
islands. By inspecting the ^Pc-plane (see the third panel in top row 
in Figs. |7|9| andpT|) we can also detect the second bifurcation of the 
zero-eccentricity equilibria which is associated with the apparence 
of a saddle at the origin accompanied by two unstable elliptic points 
close to ^1 ~ 0. These solutions can be classified as members of 
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Figure 10. Levels of the secular energy (i^ec) at the iP^-plane defined through {e\ sincoi , 62 sinc02) with coi , CO2 = ±7r/2, for a = 0.2 and n = 2.0 and for 
varied values. From the left to the right and from the top to the bottom: 0.06,0.13,0.21,0.29 (the top row), 0.37,0.45,0.53,0.61 (the middle row), and 
0.69,0.76,0.84,0.92 (the bottom row), respectively. Shaded regions illustrate mutual inclinations in prescribed ranges (lower inclination — lighter shade, 
larger inclination — darker shade). For each fixed value of J^, there are two levels of Imut which are marked in subsequent panels: (25°, 35°), (45°, 55°) 
(65°, 75°), (75°, 85°) for the upper row; (90°, 100°), (100°, 110°), (110°, 120°), (120°, 130°) for the middle row; (130°, 140°), (140°, 150°), (150°, 160°), 
(160°, 170°) for the bottom row, respectively. 
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Figure 12. A comparison of families of stationary solutions found with the octupole expansion of the secular Hamiltonian (green and blue filled circles) and 
with the help of quasi-analytic method (black and red filled circles). Parameters are for a = 0.333 and /j^2. Large filled circles are for stable solutions, and 
smaller circles are for unstable equilibria. The red and green filled circles are for equilibria beyond the formal convergence limit of i^ec expanded in a. 



family lib because they emerge in quadrant II of the 5 -plane. The 
saddles visible in ^Ps (Fig.jsJ represent members of family Ilia [see 



also Fig. 10 for the second pair of (a,/i)]. 



Thanks to the semi-numerical algorithm, we can follow the 
zero-eccentricity family not only for larger a, up to 0.667, but also 
in wide ranges of the mass ratios, between 0.25 and 2. Stability 
diagrams in Figs. pJjfTS] illustrate bifurcations of equilibria for dif- 
ferent parameter pairs. Obviously, the zero-eccentricity solutions 
appear for all combinations of the primary parameters, moreover, 
the evolution of this family is very complex. Bifurcations of family 



"0" lead to new classes of equilibria and qualitative changes of the 
i^ec topology as seen at the 5 -plane. 



4.2 Family la 

Family la appears for a limited range of (a,/i) in quadrant I of the 
5-plane. We detected it for ^ ^ 1 (more massive inner planet). This 
type of stationary solutions is characterized by small ei and a range 
of 62 between and a value permitted by the equation of the col- 
lision line. It emerges from a bifurcation of the zero-eccentricity 
solution in the range of large and is always unstable. 
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Figure 13. Families of stationary solutions obtained for a = 0.2 and the following mass ratios: n = 2.0 - the top left-hand plots, fi^ 1 .0 - the top right-hand 
plot, ;U = 0.5 - the bottom left-hand, and fi = 0.25 - the bottom right-hand plots. Large filled circle are for Lyapunov stable (or linearly stable) equiUbria, 
smaller filled circles are for unstable equilibria, red filled circles are for solutions found in regions, where the power series of i^ec in oc would diverge. The 
stationary solutions are classified according with the quadrant of the 5-plane, in which they appear, hence the columns in each sub-group of diagrams [for 
fixed (a,^) written in the legend] are for the following (0)1,0)2) -pairs: (7r/2,7r/2) - the left column, (7r/2, — 7r/2) - the middle column, and (0,0) - the right 
column. Each sub-group of stability diagrams has panels for the mutual inclination (the top row), and for the eccentricities (the middle and the bottom rows, 
respectively). The x-axis of each diagram is labeled by and Iq. Particular families of solutions which are identified in this work are labeled with Roman 
numbers and appropriate Latin letters. See the text for more details. 
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Figure 14. Families of stationary solutions obtained for a = 0.333 and the following mass ratios: ^ = 2.0 - the top left-hand plots, ;U = 1 .0 - the top right-hand 
plot, ;U = 0.5 - the bottom left-hand, and ja = 0.25 - the bottom right-hand plots. Large filled circle are for Lyapunov stable (or linearly stable) equiUbria, 
smaller filled circles are for unstable equilibria, and red filled circles are for solutions found in regions, where the power series of i^ec in oc would diverge. 
The stationary solutions are classified according with the quadrant of the 5-plane, in which they appear, hence the columns in each sub-group of diagrams [for 
fixed {oL,jn) written in the legend] are for the following (0)1,0)2) -pairs: (7r/2,7r/2) - the left column, (7r/2, —n/2) - the middle column, and (0,0) - the right 
column. Each sub-group of stability diagrams has panels for the mutual inclination (the top row), and for the eccentricities (the middle and the bottom rows, 
respectively). The x-axis of each diagram is labeled by and io. Particular families of solutions which are identified in this work are labeled with Roman 
numbers and appropriate Latin letters. See the text for more details. 



© 0000 RAS, IVINRAS 000, 000-000 



Equilibria in non-coplanar secular problem 15 



coi=7u/2, co2=:7i;/2 



COi=7U/2, C02=-7l/2 



COi=0, 0)2=0 



coi=7u/2, co2=:7i;/2 



COi=7U/2, C02=-7l/2 



90 



io[deg] 90 



45 



1 


1:: 1 1 








IVc 












Vb- 












* f h« 1 


i 







/o[deg] 90 



io [deg] 



1:: 1 


: 1 1 


M 






^ Ila 








1 r**-i 




FT 


; 1 1 1 

Ille 












I IVb+\ 

: 1 lii^ 


L 



FT 


1 












1 1 


lib 






;; ( 


1 1 


1 


L 



0.5 



A 



0.5 



0)i=7i;/2, o)2=7r/2 



0)i=7r/2, o)2=-7i/2 



0)^=0, 0)2=0 



90 180 [deg] 90 180 [deg] 90 180 [deg] 



90 
45 



1 1 ::l 1 1 
[ 1 IVb-> 


• 


r 1 


1 



1 1 


1 1 1 






Ilia 
1 1 


IIIC+ 

1 



1 1 


;l 1 1 






illb 


TTn 




; 1 ^i^*^ 




1 








a=0.667 


• 
• 

• _ 




|i=0.5 


llld/ 










1 N 


^IIIc+ 




I Ilia 


\ 




I, 





n~n — i~r 



y mdi r 



1 1 


1 11 


Ila 




lib 

-1 1 


- 1 




1 1 


1 1 1 

/ 












;IIb 
1 1 


1 



A 



0.5 



A 



0.5 




io[deg] 90 



/o[deg] 



1 1 


1 1 


1 














^Ila 








1 1 


il 






0.5 





^2 



r 



1 :: 1 1 


1 


1 


1 








-i- 


IVa 








^^1 1 


^1 







r-m — r 



a=0.667 
1^=1.0 



Ilia 



Illf 



1 ^ 




1 


ll^Ia 






-r M 






■ lib 

: V 1 1 


^1 





0.5 



n^n — r 



iVa 
J_ 



r-m — r 



fllla 



jiif 

/H 



FT 



1 — r 



■^lla 



i^IIb 



J- 







0.5 



A 



0.5 



A 



0.5 



0)i=7i;/2, o)2=7r/2 
90 180 in [deg 



^mut 


1 1 ;;i 1 tt 


135 




90 


\ :m \ • 




45 


rT\iva 




^7 


III 1 





0)i=7r/2, o)2=-7i/2 0)^=0, 0)2=0 

90 180 /^[deg] 90 180 [deg] 




FT 


FT 

y 






lib 






Ml 1 




Figure 15. Families of stationary solutions obtained for a = 0.667 and the following mass ratios: ^ = 2.0 - the top left-hand plots, ;U = 1 .0 - the top right-hand 
plot, ;U = 0.5 - the bottom left-hand, and fi = 0.25 - the bottom right-hand plots. Large filled circle are for Lyapunov stable (or linearly stable) equiUbria, 
smaller filled circles are for unstable equilibria, red filled circles are for solutions found in regions, where the power series of i^ec in oc would diverge. The 
stationary solutions are classified according with the quadrant of the 5-plane, in which they appear, hence the columns in each sub-group of diagrams [for 
fixed (a,^) written in the legend] are for the following (0)1,0)2) -pairs: (7r/2,7r/2) - the left column, (7r/2, — 7r/2) - the middle column, and (0,0) - the right 
column. Each sub-group of stability diagrams has panels for the mutual inclination (the top row), and for the eccentricities (the middle and the bottom rows, 
respectively). The x-axis of each diagram is labeled by and Iq. Particular families of solutions which are identified in this work are labeled with Roman 
numbers and appropriate Latin letters. See the text for more details. 
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4.3 Families IVa, IVb-, IVb+, Illb — the L-K resonance 

We have already seen that the L-K equilibrium appears in quad- 
rant IV of the 5-plane (family IVa) and is tightly related to family 
"0" because it emerges from its "first" bifurcation on the J^-axis. 
The second L-K bifurcation leads to family IVb- associated with 
saddles in quadrant IV and to a pair of a saddle (IVb-i-) and an el- 
liptic point (family Illb) in quadrant III. These structures are partic- 
ularly well seen in the middle row of Fig.[TO| Stability diagrams in 
Figs. pJ}fT5] tell us that equilibria of families IVb- and IVb-i- are lin- 
early stable while solutions of family Illb are unstable. Curiously, 
equilibria of family IVb-i- might be identified with non-restricted 
case of the L-K resonance characterized by librations of coi around 
±7i/2 with possible simultaneous librations of AG5 around (or, 
equivalently, CO2 around 

To examine more closely the secular dynamics in the regime 
of the classic L-K resonance (family IVa), we computed the 
Poincare cross sections for the secular Hamiltonian having two de- 
grees of freedom. These cross-sections are obtained by integrating 
the equations of motion over a few Myr time-scale. The parame- 
ters are selected as for the HD 12661 system (see the caption to 
Fig.[T6|. The cross-section planes are chosen as follows: 

El = {ei cosACD x ei sinACO}, E2 = {ei cos2cOi x ei sin2cOi}. 

The surface of section Ei is defined by coi =n/2 (d(Oi/dt < 0), 
and the plane E2 by AG5 = 7i {d/S^/dt < 0), respectively. The top 
panels of Fig.[T6]are for the Ei -plane, bottom panels of Fig, are 
for the E2 -plane. The cross sections are computed for energy curves 
in the neighborhood of the L-K quasi- separatrix: panels in the left- 
hand column are for the initial conditions lying on the energy level 
within the quasi-separatrix curve encompassing the L-K resonance 
center, panels in the middle column are for the energy level corre- 
sponding to the quasi-separatrix curve, and the right-hand panels 
are for the energy level encompassing the quasi-separatrix. In the 
cross-sections, we can detect a few high-order secular resonances, 
invariant curves representing quasi-periodic orbits and relatively 
large regions of chaotic motions. The appearance of chaotic dy- 
namics in the 3-D problem is the new feature as compared to the 
co-planar dynamics. We recall that in the later case, the averaging 
leads to one-degree of freedom integrable system. 

Actually, the smooth invariant curves seen in the Poincare 
cross sections assure us that quasi- analytic averaging makes it pos- 
sible to derive very precise numerical solutions of the secular equa- 
tions of motion. The Poincare cross-sections, although obtained by 
complex algorithm relying on the numerical integration of the mean 
Hamiltonian and the equations of motion derived by numerical dif- 
ferentiation, make it possible to study the dynamics of the secular 
system in detail, and in the whole permitted range of the orbital 
parameters. 

4.4 Families Ila and lib 

Equilibria of family Ila are very special, because they are found 
over the collision line of planetary orbits, hence beyond the formal 
limit of convergence of the expansion of i^ec in a. We call them 
the chained stationary configurations because both secular orbits 
are connected like the links of a chain (see Fig. [TT] for an illus- 
tration). These solutions appear at small mutual inclinations /mut- 
In spite of large eccentricities, the mean orbits cannot cross each 
other thanks to a particular spatial orientation. These equilibria are 
generically stable because they correspond to the maxima of i^ec- 
That can be seen in the stability diagrams (Figs.[T3}{T5]). With in- 




Figure 17. An example stationary configuration of family Ila computed for 
a\ = 0.83 au, a2 = 2.56 au, e\ = 0.325, ^2 = 0.7204. The mutual inclination 
of orbits in this case is /mut ~ 24° . 

creasing J^, family Ila emerges close to the collision line and then 
"moves" towards the border of the 5-plane. Curiously, the eccen- 
tricity e2 of equilibria Ila spans moderate and large values, and this 
family can be detected for all pairs of parameters analyzed in this 
work. The most prominent example of family Ila is shown in sta- 
bility diagrams computed for a = 0.667 (see Fig. [15]). We stress, 
that these solutions are non-resonant. Family lib is characterized 
by very small eccentricity of the inner planet and large mutual incli- 
nation of orbits (let us recall that it it may emerge for coi = CO2 = 0). 
This family appears at bifurcational inclination and can be re- 
vealed by the quadrupole theory (Krasinsky 1974 ) . It is unstable — 
although it appears in the 5 -plane as an elliptic point. In fact, the 
secular Hamiltonian is not sign definite function in its neighbor- 
hood. The linear stability analysis reveals complex eigenvalues of 
the linearized equations. 

4.5 Families of quadrant III of the 5-plane 

For small J^, equilibria appearing in quadrant III of the 5-plane can 
be associated with bifurcations of the zero-family solutions. These 
are equilibria of family Ilia at a saddle seen in the fP^-plane (see 
two last panels in the top row of Fig. |8]) and appearing closely to 
the ^1 = axis, and they are always unstable. Other solutions in 
quadrant III are associated with a bifurcation of the L-K resonance 
in the regime of large ei (families Illb, IVb-i-). We note, that in 
our "taxonomy", the name of IVb-i- for the saddle associated with 
elliptic point Illb is justified by the fact that this point can move be- 
tween quadrant IV and quadrant III (see a sequence of panels in the 
middle row of Fig.[To]). Equilibria Ilia appear for different critical 
values of (or mutual inclination /q) than solutions of family lib. 
This can be particularly well seen in Fig.[T5]for large a. Moreover, 
the quadrupole order theory predict that these families appear for 
the same /q. 

In the regime of large J^, a plethora of quadrant III solutions 
appears. We classify them with symbols IIIc, IIIc-, IIIc-i-, Illd, Ille, 
Illf (see, e.g., the right-hand panel of Fig. [3]) and they are associ- 
ated with large mutual inclinations of orbits and e2 ^ I, when the 
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Figure 16. Poincare cross-sections (coi = 7r/2, dodi/dt < 0) in the top row, (An3 = %, dAtn/dt < 0) in the bottom row, computed for 3D configuration of two- 
planet system with orbital parameters corresponding to the best-fit parameters of the HD 12661 planetary system: mo = 1.07 M©, mi = 2.3 mj, m2 = 1.57 mj, 
ai = 0.83 au, a2 = 2.56 au. = 0.085678. The plots in columns from the left to the right are for the following secular energies: j = -5.106073 x 10~^, 
-E^^ = -5.106139 X 10"^ 'Ecj = -5.106490 x 10"^ (in canonical units of IMq, 1 au, 1 yr and k = 27r). These energies are chosen in the neighborhood of 
the Lidov-Kozai resonance. Red dots are for chaotic motions, blue dots mark quasi-periodic solutions. See the text for more details. 



energy plane become disconnected. Some of these families can be 
linearly stable (see the stability diagrams in Figs. pj}fT5] l. In gen- 
eral, the numerical continuation of these families is very difficult 
because they evolve close to the boundary of permitted motions in 
the 5-plane, and in the regime of large mutual inclinations and ec- 
centricities. Then the numerical procedure sometimes fails, due to 
not precise enough determination of the second order derivatives, 
and that may also explain some gaps in the family curves which are 
present in the stability diagrams. Also problems with the continu- 
ation of these families hinder precise and proper identification of 
some solutions. 



5 CONCLUSIONS 

The semi- analytical averaging is a powerful technique helpful to 
reduce limitations of the analytical theories. Although it relies on 
purely numerical algorithms, its solid theoretical background is vi- 
tal for the interpretation and understanding of the results of nu- 
merical experiments. In this work, we demonstrate that it makes 
it possible to study the secular dynamics of two-planet system in 
wide ranges of semi-major axes and masses ratio. The appropriate 
scaling of the problem parameters helps us to represent the phase 
space of the secular system globally. Assuming non-resonant con- 
figurations, and that orbits are distant enough from collision zones, 
the averaged system may be reduced to two degrees of freedom. 
Hence, to carry out the analysis, we can apply geometric tools, like 
the representative planes of initial conditions, Poincare cross sec- 
tions, continuations of stationary solutions with respect to parame- 



ters, which are very helpful to understand the structure of the phase 
space and, in turn, the long-term behavior of the planetary system. 

Equilibria are the basic class of solutions which can be inves- 
tigated with relatively simple tools. Our analysis reveals a number 
of families of stationary solutions in the 3D secular problem of two 
planets. To the best of our knowledge, some of them are yet un- 
known in the literature and are related to unusual orbital configura- 
tions. For instance, we found the so called chained stationary con- 
figuration which are non-resonant, can be found in the regime of 
small mutual inclinations and large eccentricities, and are located 
over geometric collision line of orbits. In spite of such extreme dy- 
namical situation, these secular solutions are Lyapunov stable and 
exist in wide ranges of semi-major and masses ratio. Simultane- 
ously, such orbital configurations prohibit application of analytical 
methods relying on power series expansion of the perturbations. 
The semi-analytic averaging helps to generalize analytical results 
obtained for low-order expansions of the secular Hamiltonian. 

We obtained some interesting results regarding the Lidov- 
Kozai equilibrium in the non-restricted problem. We found that this 
resonance may be associated with librations of AG5 around n in the 
neighboring trajectories (not only with librations of the inner peri- 
center around ±7i/2). These librations are possible for relatively 
large ratios of semi-major axes and planetary masses. We found 
that the L-K resonance may also appear in the regime of large ec- 
centricity of the inner orbit, and then it would be associated with 
librations of coi around ±n/2 with simultaneous librations of AG5 
around 0. The parametric evolution of the L-K resonance is related 
to the stability of the zero-eccentricity equilibria. There is a link be- 
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tween bifurcations of this family (and changes of its stabiUty) with 
an appearance of new famihes of stationary solutions in other parts 
of the phase space (hence, with changes of its global topology). 

Our work illustrates qualitatively different view of the 3D dy- 
namics as compared to the coplanar configurations. It is already 
known ( [Michtchenko & K^ lhotra 2004) that the non-resonant, 
coplanar systems of two point-mass planets fulfilling the averaging 
theorem and interacting through Newtonian forces are integrable. 
The secular dynamics of such systems are basically trivial and may 
be reduced to one degree of freedom. Under the same assumptions, 
the spatial configurations may exhibit strong chaos and extremely 
complex secular phenomena. 

In the approximation of Newtonian, point-mass interactions, 
the dynamics depend on masses and semi-major axes ratios, hence 
our results are valid both for planetary systems with small planets, 
as well as for systems comprising of brown dwarfs or even sub- 
stellar companions. However, the dynamics of real systems may 
strongly depend on the magnitude of the mutual interactions. More- 
over, many stable equilibria are found for large values of J^. Ac- 
cording to the notion of AMD (Laskar "2000 ), in such cases the 
real configurations are unlikely long-term stable even close to sec- 
ularly stable equilibria. In that sense, the results of stability anal- 
ysis may be too optimistic. We skip a study of such effects in this 
work, because it would make the paper necessarily very lengthly. 
We should keep in mind that introduction of relativistic, tidal, and 
stellar quadrupole-moment perturbations, may affect the secular 
dynamics dramatically ( [Migaszewski & Gozdziewski|2009b|a] ). 

The quasi-global technique applied in this paper has been 
proved an efficient and effective tool for the analysis of the sec- 
ular 3-D model. Nevertheless, we learned from the work that the 
quasi-analytical approach is not a perfect tool. Due to limitations 
of the numerical algorithms, the continuation of families of equi- 
libria and analysis of their stability is particularly difficult when 
we reach limits of permitted motion or i^ec is a weakly varying 
function. We can also overlook some solutions. Unfortunately, that 
leave us sometimes with open questions. 
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